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Abstract 

Asymptotic expansion of the second-order linear ordinary differential equation + k 2 f(z)^f — 0, in which 
the real constant k is large and / = 0(1), can be carried out in the manner of Liouville and Green provided 
/ does not vanish. If / does vanish, however, at zq say, then Liouville- Green expansions can be carried 
out either side of the turning point z = zq, but it is then necessary to ascertain how to connect them. This 
was first accomplished by Jeffreys, by a comparison of the differential equation with Airy's equation. Soon 
afterwards, the situation was found to arise in quantum mechanics, and was discussed by Brillouin, Wcntzel 
and Kramers, after whom the method was initially named. It arises throughout classical physics too, and 
is encountered frequently when studying waves propagating in stars. This brief introduction is aimed at 
clarifying the principles behind the method, and is illustrated by considering the resonant acoustic-gravity 
oscillations (normal modes) of a spherical star. 
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1. Introduction 

It is common, and often apparently most straightfor- 
ward, to approach problems in macroscopic stellar physics 
by solving the governing equations numerically. However, 
simple analytical techniques are often more revealing. For 
the latter it is usually necessary to idealize the situation 
in hand, sometimes grossly so, to render the equations 
tractable. The outcome can reveal the phenomenon un- 
der study in a very different light from that provided by 
specific numerical examples. In particular, because in 
some respects analytical results are more general, appreci- 
ation of their structure can guide one more easily towards 
greater understanding. Even viewed merely as a diag- 
nostic tool, that understanding, and sometimes just the 
analytical solutions themselves, have proved to be useful 
in the past simply for finding errors in numerical com- 
putations, both before and after publication. Numerical 
and analytical investigations are complementary, and each 
would be much poorer without the other. 

This is the first of a short series of invited articles in- 
tended to illucidatc some of the analytical techniques that 
are used in macroscopic stellar physics. It discusses one 
of the most useful techniques for studying the wave-like 
solutions of ordinary linear differential equations of sec- 
ond order: namely, the so-called Liouville-Green expan- 
sion combined with the method of Jeffreys for connecting 
solutions across turning points, more commonly known 
as the WKB approximation. The method is presented 
and, I hope, made plausible, giving a slight flavour of the 
background arguments without resorting to mathematical 
proof, with the aim of guiding those not already conver- 



sant with the method towards its proper use. 

2. The Liouville-Green expansion 

In modern times, the term WKB approximation has 
commonly been assigned to what is more properly re- 
garded as the leading term in the Carlini-Liouville-Green 
(LG) expansion of solutions of the second-order ordinary 
linear differential equation 



d 2 * 



+ k 2 f(z)V = Q 



(1) 



in which k is real (positive) and large, and / = 0(1) is 
slowly varying, in the sense that fc>|J?J 1 | = |-dln//dz| . 
Such equations arise throughout classical physics in de- 
scribing waves, and asymptotic approximations to their 
solutions appear to have been considered first by Carlini 
(1817), actually in a study of a problem in celestial 
mechanics, and subsequently by Green (1837) and by 
Liouville (1837). 

One of the simplest physical examples is that of small- 
amplitude adiabatic acoustic waves in an otherwise homo- 
baric fluid, whose linearized governing equations are 

P0q± = ~W, (2) 

p' = Sp = c 2 Sp = c 2 (p' + I.Vpo) = -poc 2 div£ , (3) 

where p = po + p' is pressure, p = po + p' is density and £ 
is displacement, the zero denoting equilibrium value, the 
prime, Eulerian perturbation, and S denoting Lagrangian 
perturbation. The square of the sound speed in the un- 
perturbed fluid is given by c 2 = 7iPo/po> where 71 is the 
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first adiabatic exponent. The unperturbed pressure po is 
constant, and in this simple case po and c 2 are considered 
to vary with only a single Cartesian coordinate z, and 
not with time, t. The divergence of equation (2) may be 
combined with equations (2) and (3) to yield the equation 



1 d 2 5p 
c 2 dt 2 



= 



(4) 



where H p 1 n = ~ Vlnpo defines the density scalcheight H p , 
and n is a unit vector in the z direction. In the special case 
in which the disturbance Sp varies only in the z direction, 
equation (4) reduces to the Klein-Gordon equation 



d 2 X 
dt 2 

where \ : 



2 2 

w c X = c 



d 2 x 
dz 2 

Sp and 



AH 2 



dH„ 
1-2 p - 

dz 



(5) 



(6) 



Because po and c 2 are constant in time and are inde- 
pendent of the co-ordinates perpendicular to n, equation 
(4) admits oscillatory solutions with frequency u) of the 



form 5p(x,t) = 3? p^{z) 
real part and \I> satisfies 



d 2 * 

dz 2 



i(kj 



* = 



ji ) 



where 5ft denotes 



(7) 



in which fcj_ = |fcj_|, where k± is a constant wavenumbcr 
component perpendicular to n. Equation (7) is essen- 
tially of the form of equation (1), provided that lo is high 
enough. 

If po and c were actually constant, lo c would vanish and 
equation (7) would admit solutions ^ = \I/oe ±l,C||Z , where 
is a constant amplitude and k\\ = \Jlo 2 /c 2 — k^_. This 
represents a wave travelling with speed := to/ \ k\ = c 
in the direction k = (fej_,±fc||), provided that fc|| is real, 

because the phase <fi = k.x — Lot is invariant in a frame 
of reference moving in the direction of k with speed v^, 
namely with x = xq + v^tk, where k = \ k\ 1 k. The quan- 
tity i>0 is called the phase speed. 

This solution motivates the LG expansion. Suppose 
first that / > in equation (1). If c 2 and lo 2 in equation 
(7) vary 

slowly with z one can write (lo 2 — lo 2 ) /c 2 — k 2 _ = k 2 K 2 (z), 
where k is constant and k = f 1 / 2 is a function whose mag- 
nitude is of order unity and whose scalcheight H K is much 
greater in magnitude than k . In other words, the equi- 
librium, background, state varies only little over the char- 
acteristic lengthscale of variation of the wave. It is 
only under such circumstances that the concept of a wave 
is readily interpreted. 

Thus, one poses a wave- like solution to equation (7) by 
setting 



z<fi(z,t) 



i(kip(z)—Lut) 



(8) 



regarding k as a large parameter. It is usual then to ex- 
pand ip in inverse powers of k; here I shall afford myself 



the flexibility of expanding both ip and \I>o : 

i> = i> + k~ Vi + k ~ 2 ip2 + ■■■ , 

*o = *o,o + fc _1 *o.i + fc _2 *o,2 + ... . (9) 

Substituting expressions (8) and (9) into equation (7), 
with k 2 K 2 in place of [lo 2 — lo 2 ) /c 2 — k\, and equating to 
zero terms of like order, yields the sequence of equations: 



2# tf'+(i<-2V> V 1 )*o,o = 



(10) 
(11) 



2i^ * , o,i + Wo ~ 2^0*0,1 + *o,o + 2*^i*o,o 

+ (i<-^ 2 - 2^)*o,o = 0, (12) 

etc., 

where here and henceforth the prime denotes differentia- 
tion with respect to the argument. These equations can 
be solved successively. 

The first is equation (10), whose solution is 



± / ndz 



(13) 



One now encounters in equation (11) a nonuniqucness, 
because we have two new functions, ^o,o and ip±, with 
only one equation to determine them. This is a result of 
the flexibility I introduced by expanding the two functions 
and tp in a representation of only a single function 
iff, so there is redundancy. One is therefore at liberty 
to choose either one of those functions, or any relation 
between them, as one wishes. How one does that can be 
regarded as a matter of convenience, and first I make the 
common choice of setting ip\ = and solving the resulting 
equation for the leading-order amplitude, ^0,0 ; yielding 



#0,0 OC^' V2 <X K- 1 ! 2 



(14) 

the signs of proportionality indicating that one can multi- 
ply the functions by any constant. One can continue, but 
evidently the equations are beginning to become some- 
what cumbersome; rather than trying to write down a gen- 
eral procedure, it is preferable to tailor one's way through 
the complexity, using prudence to discern the way. 

Of course one could cut the cackle by setting ip n — 0, 
n > 0, yielding 2tP'% n + ip"^ Q , n - i*o,n-i = 0, V = V>o, 
whence 



#0,n = pK 



(15) 



this has the formal advantage of generating the entire se- 
quence of functions in one compact formula. However, 
the sequence does not always converge. Evidently, if one 
expands only $q and if ip does not describe the phase ade- 
quately, the (real part of the) exponential e %k ^ will vanish 
in the wrong place, and although "J might be forced to 
have a zero in the right place in an attempt to rectify 
that, it cannot remove a zero in 5ft(e lfc ^) that is in the 
wrong place without itself being singular at that point. 
Therefore, the radius of convergence of the expansion is 
likely to be smaller than it need be. Alternatively, one 
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could expand only ip, as was common in the early days. 
The procedure is algebraically more complicated, but it 
tends to deliver a more robust approximation. Alternate 
terms are real and imaginary, the latter accounting for the 
variation of the amplitude. In this respect, it is interesting 
to observe that if instead of setting ipi = in equation (11) 
one sets ^0,0 = L then ip[ = iipQ /2ip' ; whence tpi = |ln«, 
and 



±?fc / K,dz— ilnK —1 ±ik I ndz 



(16) 



which reproduces the solutions implied by equations (13) 
and (14). Proceeding further, one obtains the second- 
order correction to the phase, which is given by equation 
(17) below; then the third-order term provides the second- 
order correction to the amplitude, which is identical to 
equation (18). 

The algebraic technicalities are eased by expanding 
both ip and "Jo- One encounters nonuniqueness at each 
stage in the sequence of equations (10)-(12), and it is of- 
ten expedient to alternate between expansions in ^>q and 
ip, which is not surprising given that the terms in the 
phase expansion alternate between being real and imagi- 
nary when ^0 is held fixed. Thus, one can set ^0,1 = in 
equation (12) and obtain 



V; 2 = ±i/^ 1/2 



dt 2 ^ ) 



(17) 



and then set ips = in the subsequent equation in the 
sequence to yield 



H 2 

iTr 1 ..-2 u „-l/2 

*o, 2 = - 3 « cL? 



(18) 



and so on. Equations (17) and (18) are somewhat sim- 
pler to derive by this procedure than are their counter- 
parts in the pure phase expansion. The expressions for 
ip2 and ^0,2 are essentially first and second derivatives of 
K -1 / 2 respectively; higher-order 'corrections' to the solu- 
tion contain yet higher derivatives, which gencrically au- 
gur eventual divergence: the expansion is asymptotic and 
must be terminated at some order. 

If one needs to develop connexion formulae of higher- 
order than those presented in §5 it is prudent to arrange 
for ^ ,n and ip n +i to be related in such a way as to ensure 
that the Wronskian , J , + , I''_ — \K of the approximations 
and to linearly independent solutions of equation 
(7) is constant (Froman and Froman, 1996), as it is for the 
exact solutions, and as it is also for each of the asymptotic 
representations (16) and (19). 

Terminating the LG expansion at equation (14) pro- 
vides the simplest, leading-order, representation. It re- 
sults from solving just the first pair of equations in the 
sequence beginning with equations (10)— (12) . It is that 
formula which nowadays is often called the WKB approx- 
imation. 

In cases where / < in equation (1), k 2 < in equation 
(7), solutions can be obtained by exactly the same proce- 
dure as the wave-like solutions, provided \f(z)\ remains of 
order unity so that the ordering of the sequence of equa- 
tions (10)-(12) etc. is preserved. The outcome in the 



so-called WKB approximation is again given by equation 
(16), which is better rewritten: 



| K |-l/2 e ±kJ>|d, 



(19) 

The solutions therefore vary exponentially. Such 'waves' 
are termed evanescent; evidently they do not propagate 
in the z direction. 

Near a point where k 2 vanishes, the ordering of the 
sequence of equations (10)-(12) etc. is not preserved, and 
the approximation cannot be used. 

3. Normal form 

Equation (1) is in what is called normal form, having no 
term in which "J is singly differentiated. It might appear 
at first sight that to add such a singly differentiated term 
would be of no serious consequence, because the formal 
LG expansion could still be applied to the more general 
equation. However, to do so is not prudent, for unless 
one is very careful indeed, and perhaps even if one is very 
careful, one risks ending up with a representation of the 
solution whose domain of applicability is more restricted 
than it need necessarily be. To illustrate the point let us 
consider again an equation with constant coefficients: 

d 2 * 



dz 2 



2?7^ + fc 2 * = 
dz 



(20) 



with 7] = 0(1). Its solutions are exp(— r\z ± i\/k 2 — rj 2 z). 
Applying the LG expansion in the usual way, up to 0(k), 
yields the WKB equations: 

ip'o = ±k , 



(21) 
(22) 

(23) 

This approximate solution is valid only for a range d of z 
satisfying d^2k/r] 2 . 

If, on the other hand, one writes *B = e~ nz u(z), then 



*o,o 
whence 



= -v 



-r}Z±ikz 



— + (k 2 - V 2 )u = 0, 



(24) 



which is in normal form. Now the WKB-approximate so- 
lution is actually exact. Of course, one could continue 
the direct LG expansion of equation (20) to higher order, 
which simply provides the expansion of \/k 2 — r\ 2 z in in- 
verse powers of k 2 , but it takes the (formal) solution of 
two of the differential equations in the analogue of the se- 
quence (10)-(12) etc. for each term, which is a lot of work. 
Therefore it is undoubtedly at least expedient first to cast 
equations with singly differentiated dependent variables 
into normal form before proceeding with the expansion. 

Of course, it could be said that the demon- 
stration with an equation with constant co- 
efficients proves nothing about expansions 
of equations with non-constant coefficients. 
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That is not wholly the case, however, because pro- 
vided the coefficients vary smoothly, so do the solutions, 
and one can say that the casting into normal form is 
advantageous at least for equations whose coefficients are 
close to being constant. Furthermore, we have a great 
deal of experience with equations of this kind, and we 
know that under a wide variety of circumstances these 
asymptotic techniques based on small departures from 
constancy work much better when the departures are not 
so small than perhaps one might feel they ought. This 
is not mathematical proof, but pragmatism, based upon 
which I strongly recommend the taking of the trouble 
to formulate the problem sensibly at the outset. To do 
so is unlikely to cause a deterioration in the eventual 
outcome, and I baldly assert that it is actually very likely 
to provide substantial improvement. 

Second-order linear ordinary differential equations can 
always be cast into normal form, just as they can always 
be case into self-adjoint form. If one has, for example, the 
equation 



+2?7— + fc K y = 
dz z dz 



(25) 



where k is constant, and 77 and k are functions of z, and if 
one wishes to preserve the independent variable 2, then 
one writes y = substitutes into the differential equa- 
tion, and simply chooses u in such a way as to make the 
coefficient of /dz vanish. The result is u — exp(— J rjdz) 
and 



(26) 



It is also sometimes desirable to transform the inde- 
pendent variable into something more natural, such as 
acoustic radius (namely, sound travel time from the cen- 
tre of the star), for example, if one were studying stellar 
acoustic waves. That alone would destroy the nor- 
mal form, if the equation had been in that form in the 
first place. However, one can again start from equation 
(25), but first express it in terms of a new independent 
variable x, and once again write y = w^, this time with 
u = (dz/dx) 1 / 2 exp(— J rj^dx), to obtain 



d 2 fr 

da; 2 



where 



v = w 



v(x) 



* = , 



dw 
dx 



w(x) = 1] 



dz 
dx 



Id, f dz 
In — 

2 dx \dx 



(27) 



(28) 



In the case of waves propagating in the z direction, it can 
be useful to replace z by the mass variable q = J p dz. 
Waves in either a homobaric fluid or in a plane-parallel 
atmosphere stratified under constant gravity g satisfy 



d 2 Sp _ ^d 2 5p 



dt 2 



dq 2 



(29) 



where c = pc, which establishes a transformation between 
the Klein-Gordon equation (5) and the wave equation 
(29). By applying the WKB approximation to the usual 



wave-equation analysis it is straightforward to demon- 
strate that an arbitrary infinitesimal disturbance dp prop- 
agates at a rate c along the q coordinate in cither direction 
without significant change of shape and with amplitude 
varying as c 1 / 2 , provided that the scale of variation of 
that disturbance is much less than the scalchcight of c. 
If the disturbance varies sinusoidally with time, the wave 
equation (29) reduces directly to the form (1). 

4. Critical acoustic cutoff frequencies 

Before proceeding to a discussion of the JWKB approx- 
imation, which I intend to illustrate with the problem 
of determining asymptotic properties of acoustic-gravity 
waves, I pause for a moment to discuss the so-called crit- 
ical cutoff frequencies. They represent the frequency 
beneath which 

a wave cannot propagate, although, when described 
in 

these physical terms, it must be appreciated that they are 
not uniquely defined. That is simply because the concept 
of propagation itself is not well defined. As I alluded at 
the outset, even the idea of a wave is itself an asymptotic 
concept, and is not easily interpreted unless the scale of 
variation of the background state is substantially greater 
than the magnitude of the inverse wave number - the very 
conditions under which the approximations described here 
are designed to be used. It has been a consequence of 
the resulting imprecision that some apparently unsuspect- 
ing workers have been too careless in writing down what 
should have been precise equations to describe a physi- 
cal situation, naturally under idealized, yet well defined, 
conditions, and have so degraded their inferences unnec- 
essarily. 

The critical frequency that is perhaps the most famil- 
iar in physics is the plasma frequency, lu p . Indeed, its 
very existence is the reason why plasma is so named (I 
refer here to an ionized gas, not to the liquid component 
of blood, although the basis for the latter's appellation 
is essentially the same.) The advantage of this example 
over that of acoustic-gravity waves is that it can be con- 
sidered for an infinite uniform plasma (in the absence of 
an imposed magnetic field), in which Langmuir waves of 
frequency w satisfy the equation 



to 



d 2 ^ | 

dz 2 c 

„2 OT1 j ,.,2 _ „ Jl 



(30) 



with c and u> 2 = n e e /m e eo each being constant. In that 
case the critical cutoff frequency lo p is quite well defined 
in the physical terms I used above, as is the concept of 
propagation: high-frequency waves, with u> > u> pi prop- 
agate in the z direction (which, of course, is arbitrary) 

with wavenumber k = yjuj 2 ~~uj 2 /c, whereas temporally 

periodic disturbances with u) < uj p are evanescent, having 



an e-folding length c/y uj 2 — u 2 . 

I might mention, in passing, that Langmuir waves have 
been likened to the Jeans waves in an essentially infi- 
nite uniform self-gravitating fluid, a situation with which 
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many astronomers are more familiar. These waves are 
said to satisfy 



d 2 ^ 

dz 2 



-* = 



(31) 



(they actually do, but only approximately, for to avoid 
the so-called Jeans swindle the background state must 
vary, slowly, in either space or time, as Jeans himself rec- 
ognized), where w 2 = AitGpo. However, with respect to 
scales much smaller than the Jeans length c/wj, the grav- 
itational term is negligible, and equation (31) loses its 
cutoff. It is only when po or c 2 varies with z that a cut- 
off is reintroduced, but then the spatial variation tends to 
cloud the idea of propagation. The resulting imprecision 
is emphasized by comparing the Klein-Gordon equation 
(5), which contains an explicit cutoff frequency u> c , with 
the precisely equivalent wave equation (29), which does 
not. The issue at stake, if one wishes to retain the phys- 
ical picture, is to define, in a spatially varying medium, 
where the wave can propagate and where it cannot. I hope 
it is evident now that that is unlikely to lead to a univer- 
sally well defined end. But it can be well defined under 
restricted circumstances, restricted in a sense that I shall 
explain below. And provided that one confines oneself 
consistently to the circumstances in which one has chosen 
to pose the problem, a workable definition can emerge, 
and with it the chance of a correct solution to the prob- 
lem in hand. Indeed, it is to the task of determining how 
the solutions of equation (1) in this uncertain hinterland 
connect the well determined representations (16) and (19) 
that the JWKB procedure is addressed. 

Before proceeding, permit me to digress on the issue 
of the choice of dependent variable. It is not uncom- 
mon to work with the component £|| := n.£ =: £ z of the 
displacement, or, cquivalently, velocity, in the direction 
n of variation of the background state, rather than the 
(Lagrangian) pressure perturbation. In the very simple 
case of the pure acoustic wave considered in §2, one first 
separates the parallel component £|| from the component 
of £ perpendicular to n, and then eliminates the perpen- 
dicular component and Sp from equations (2) and (3). The 
procedure is straightforward, and yields 



dt 2 



(32) 



in which V j_ is the V 2 operator in the plane perpendicular 
to n;also H ll = (— n.Vhryi) -1 and H c 2 = (— n.Vlnc 2 ) 1 
are the scaleheights of 71 and c 2 , respectively. This equa- 
tion is rather more complicated than equation (4). In 
particular, it is of higher order, although after effecting 
the separation 



£l! = 3? *(z)e 



i(kj 



.x- 



i)t) 



the resulting ordinary differen- 
tial equation satisfied by &(z) is of only second order, as 
is equation (7). The corresponding equation pertaining 
to a background state that has non-Cartesian symmetry, 



such as spherical symmetry, is yet more complicated. It 
provides justification for working with an intrinsic scalar, 
such as Sp, rather than the component of a vector. I re- 
mark also that in this simple case in which p$ is constant, 
the Eulerian and Lagrangian pressure fluctuations are nu- 
merically the same; I work formally with Sp rather than 
p' because that generalizes more easily to the situation in 
which the equilibrium state is stratified under gravity. 

Let us consider now the adiabatic acoustic-gravity 
waves in a spherical star. The governing equations in 
the Cowling approximation (an approximation obtained 
by neglecting the Eulerian perturbation to the gravita- 
tional potential) are 



dr 



L 2 g 



a- 



dSp 
dr 



9 



gpF 



1 



L c 



8p 
pc 2 







(33) 



(34) 



(e.g. Cough, 1993) with respect to spherical polar co- 
ordinates (r,9,<p), where £ is the vertical component of 
displacement from which has been factored a spherical 
harmonic function Y™(9,4>), Sp continues to be the as- 
sociated Lagrangian pressure perturbation, but now with 
Y t m factored out, g is the local acceleration due to gravity, 
with scaleheight H g , and 



F = 



9 



r L 2 q 
2+ u ^ 

Hg LOT 



(35) 



which I call the f-mode discriminant; L 2 = 1(1 + 1), where 
/ is the degree of the spherical harmonic that describes 
the angular variation of the cigenfunctions. The other 
variables continue to retain the meanings I assigned to 
them earlier. One can now eliminate £ from equations 
(33) and (34) to yield a single second-order differential 
equation for Sp, which, after reduction to normal form, 
is approximated by 



d 2 * 

dr 2 



where = (r 3 / gpF) x l 2 Sp and 
u; 2 - lo 2 c L 2 / N 2 
c 2 r 2 V <->i 2 



K 2 = 



in which 



1 



(36) 



(37) 



(38) 



is the square of the buoyancy frequency. 

Equations (36)-(38) generalize equation (7). They are 
not the exact equations to result from equations (33)-(35), 
but are what I call the planar approximation to them, 
valid as Kr — > 00. They do not include the local effect of 
spherical geometry, the only sign of sphericity that sur- 
vives being the globally geometrical representation L/r of 
the horizontal wavenumber. Including all the geometrical 
terms is straightforward (Gough, 1993); they merely add 
a little complexity to the formulae without changing those 
aspects of the mathematical structure of the equation that 
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concern us here, so I have omitted them for clarity. The 
quantity oj c , which is defined by equation (6), is what 
is called the acoustic cutoff frequency. It is not exactly a 
general cutoff frequency for \1/ in the sense that I described 
cutoff in connexion with equations (30) and (31), but is 
instead what that frequency would be for spherically sym- 
metric (L = 0) waves, uninfluenced by buoyancy. 

Equation (36) is similar to equation (1), with k = cu/c, 
c being a characteristic value of c (such as the value at 
the turning point, where / = 0), and with / depending on 
a parameter a = N/w, where N is a characteristic value 
of N. That view is appropriate for discussing acoustic 
waves. Strictly speaking, equation (36) is not precisely of 
the form (1), because a depends on fc, but evidently / 
becomes only very weakly dependent on k as k — * oo and 
the validity of the asymptotic arguments is unaffected. 
For gravity waves one takes k = L 2 N/uj and a = u>/c. 

Forgive me for emphasizing at this stage what should 
be perfectly obvious: there is no direct physical relation 
whatever between the acoustic cutoff frequency oj c and 
the buoyancy frequency N . The acoustic cutoff, which 
appears to have been discussed first by Lamb (1909), 
arises when acoustic waves cannot propagate vertically 
because the inverse wavenumber is comparable with the 
density scaleheight; consequently there is inadequate in- 
ertia on the low-density side of a compression to resist 
the inevitable acceleration of matter, thereby annulling 
too much of the pressure gradient to permit adequate 
subsequent compression of the surroundings, essential for 
causing the perturbation to propagate in a wave-like man- 
ner. The dynamics operates on the vertical component 
of the motion, and is most effective for motion that is 
purely vertical: that motion has no horizontal varia- 
tion. Buoyancy, on the other hand, exists only when there 
is horizontal variation (cf. Reye, 1872) and therefore L 
must be nonzero, as is evinced by equation (37). One way 
of regarding it is to observe that the force of gravity acting 
on a horizontally varying Eulerian density perturbation is 
not in hydrostatic equilibrium, and the unbalanced pres- 
sure gradients that are so engendered cause any typical 
fluid clement to be accelerated. That describes the pre- 
dominant dynamics of gravity waves. Confusion in the 
scientific literature between the two totally different pro- 
cesses characterized by lu c and TV appears to have arisen 
because, at least in an isothermal atmosphere with con- 
stant 71 , the formulae for the two quantities can be made 
to look somewhat similar, and, if 71 = 5/3, their values 
are almost the same, differing by only 4 per cent. That is 
no case for hiding the stark distinction between them. 

Permit me also to make another point which I hope by 
now is also quite obvious. There is a clear procedure for 
eliminating £ from equations (33) and (34) — one differ- 
entiates equation (34) to obtain an equation for d 2 Sp/dr 2 
in terms of d£/dr and £, and also dSp/dr and Sp, then 
substitutes for d£/dr using equation (33), and then for 
£ using equation (34), leaving a second-order differential 
equation for Sp - and a well defined procedure for casting 
the resulting equation into normal form, which I described 
in §2, yielding a unique dependent variable \& (to within 



an inconsequential multiplicative constant) and a unique 
equation (36). Therefore the structure of equations (36)- 
(38) is well defined, and so therefore is the acoustic cutoff 
frequency lo c , and one is not at liberty to change it. I has- 
ten to add, however, that this conclusion holds only within 
the restrictions I have imposed upon myself: namely to use 
Lagrangian pressure perturbation (or, more precisely, the 
appropriate multiple of it) as my dependent variable, and 
radius r as my independent variable. 

The full critical cutoff frequencies lo± associated with 
equation (36) for \& are easily determined by factorizing 
K 2 (Dcubncr and Gough, 1984): 



c 2 K 2 



where 



(u, 2 - W 2 )(o, 2 -u, 2 ), 



-X s 



?)± 



1/2 



and where 

<7 - LC 
r 



(39) 



(40) 



(41) 



which is sometimes called the Lamb frequency. The situ- 
ation is thus rather more complicated than it is for pure 
acoustic waves, in which buoyancy plays no part, and for 
pure gravity waves (in an incompressible fluid, in which 
sound plays no part), which I have not discussed explic- 
itly here. Nevertheless, it is apparent from equation (39) 
that solutions resembling propagating waves of the form 
(8) can be found with ip real for waves with oj 2 > and 
for waves with oj 2 < ui 2 _ , and indeed can be approximated 
by the WKB solutions (16) provided K 2 is large. If w 2 lies 
between w 2 and w 2 , then K 2 < 0, the waves can be re- 
garded as being evanescent, and can be approximated by 
the solutions (19). What is most commonly encountered 
in practice is a spatially varying yet temporally invariant 
background state, in which the frequency w of a wave is 
a conserved quantity, a property which here I take for 
granted. There the wave can encounter regions in which, 
for given lj, K 2 > 0, and therefore \&"/\l> < - the hall- 
mark of a typical wave - and regions in which K 2 < and 
vf/vp > 0. They are separated by well defined points at 
which K 2 = 0, and \S r " = 0. It is therefore convenient to de- 
fine the former regions, quite precisely, as regions of prop- 
agation (often abbreviated as propagating regions, even 
though the regions themselves do not propagate) , and the 
latter as regions of evanescence (or evanescent regions). 
They are separated by the points at which \f r " = 0, where 
the waves turn from one form to the other; these points are 
called turning points. Generally, waves in the propagat- 
ing region cannot penetrate significantly beyond a turning 
point. Yet the waves are very smooth there = 0), so 
in reality unaccounted dissipation processes cannot be in- 
voked to destroy them, in contrast to the dynamics in the 
vicinity of a critical layer (e.g. Booker and Brcthcrton, 
1967), for example, where waves are absorbed. Therefore 
isolated turning points must be points of total reflection. 
If there are two closely spaced turning points enclosing an 
evanescent region, then, of course, barrier penetration can 



No. ] 

occur, and reflection is not total. 

Although the turning points of equation (36) are defined 
precisely, they do depend on the choice of both depen- 
dent and independent variables. A different dependent 
variable, such as the displacement, £, and, more perti- 
nently, its associated counterpart S that satisfies the nor- 
mal form of the governing equation, is not in phase with 
the Lagrangian pressure perturbation, and its points of 
inflexion (together with the corresponding acoustic cut- 
off frequency) must therefore be different. Indeed, so too 
do the local vertical wavenumbers differ. But they are all 
well defined. Carrying out the procedure corresponding to 
the derivation of equation (36) yields, again in the planar 
approximation, 




the mathematical structure of which is superficially simi- 
lar to that of equations (36)-(38) and (6). The scale -H(e) 
is defined according to 

it depends on u>, rendering this formulation of the problem 
actually rather more complicated than that in terms of 
dp. However, the formula for the corresponding acoustic 
cutoff frequency ws c , defined in the sense of being the 
cutoff frequency for propagation of waves with L = 0, as 
is the cutoff frequency lo c defined by equation (6), is not 
dissimilar (aside from a sign) to equation (6), with H<a) 
being instead the scaleheight of c. I have had to adorn the 
vertical wavenumber K (and the acoustic cutoff frequency 
u> c and the buoyancy frequency N) with the subscript H 
to distinguish them from their counterparts (37) (and (6) 
and (38)), to which, for consistency, I should attach the 
subscript "P. 

A distortion of the independent variable also changes 
the locations of the points of inflexion of and H, by 
an amount which is defined by equations (27) and (28), 
although, of course, they must always occur on the evanes- 
cent sides of the locations of the nearest maxima to the 
evanescent regions. Transforming the independent vari- 
able in equations (36) and (42) to acoustic radius r = 
Jc~ 1 dr, for example, yields 

^ + (c?Xf-<)c** = (46) 
and 

^ + (^1-^5 = 0. (47) 
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In each case the square of the appropriate acoustic cutoff 
frequency is augmented by 

^-sjHfr)' (48) 

where T c — (— dlnc/dT) -1 is perhaps properly called the 
sound-speed scaletime. Equation (47) is analogous to 
a similar equation presented by Christensen-Dalsgaard, 
Cooper and Gough (1983) describing spherically symmet- 
rical (radial) adiabatic pulsations of a star with the pertur- 
bation to the gravitational perturbation included (which 
for radial waves can be cast as a second-order differen- 
tial equation) . If that equation is reduced to the Cowling 
approximation, it agrees with equation (47) with L = 0, 
as indeed it must. These equations also reduce essentially 
to corresponding forms presented by Schmitz and Fleck 
(1998) in the case when 71 is assumed to be constant. 

5. The JWKB approximation 

Consider a wave, given approximately by equation (16), 
propagating to a turning point, beyond which it is evanes- 
cent. The wave is therefore reflected, and travels back 
into the region of propagation. An interesting and im- 
portant question is: What is the phase change, if any, on 
reflection? Viewed mathematically, one has the two solu- 
tions (16) well inside the region of propagation, z < zo say, 
where zq is the location of the turning point, the positive 
and negative signs in the exponent representing the inci- 
dent and the reflected wave respectively. In the evanes- 
cent region well beyond the turning point, equation (19) 
holds; here one must choose the negative sign, because 
there is no disturbance far beyond the point of reflection. 
The question can therefore be restated thus: what com- 
bination of the two solutions (16) in z < zq match onto 
the decaying solution (19) in z > zq? Readers not inter- 
ested in the mathematical background to the answer to 
this question could skip the following subsection, save to 
accept the approximations (55)-(58) to the two standard 
solutions Ai and Bi of Airy's equation (49). 

5.1. Airy's equation 

The question posed at the end of the previous para- 
graph, in a somewhat different guise, had occupied the 
minds of mathematicians such as Stieltjes and Stokes in 
the mid 19th century. In particular, the young Stokes had 
wondered why it is that one cannot simply analytically 
continue representations (19) - actually, viewed explicitly 
as a representation of the Airy integral (50) rather than 
a representation of the solutions of the differential equa- 
tion (49) - across the turning point (to be more precise, 
around it, in the complex-^ plane, to avoid the singular- 
ity in the representation) into the representations (16). 
Although the WKB solutions (the terminology I use here 
is evidently modern, for neither W, K nor B had yet been 
born) arc invalid near the turning point, and therefore 
cannot be properly connected on the real-z axis, perhaps 
one might connect them elsewhere in the complex plane 
along some contour C chosen to be far enough from the 
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Fig. 1. The complex- 1 plane divided into three regions by 
the straight lines S~ , So and 5+ , each subtending an angle 
of magnitude 2n/3 with the others. The Airy integrals are 
defined along the infinite contours Co and C+, each of 
which asymptote to <S_ , So or 5+ at infinity. 

turning point that \kf \ ^> 1 everywhere on C. That it 
cannot be done to produce the correct solution to the 
differential equation was subsequently named the Stokes 
phenomenon. What Stokes wanted was to understand this 
matter, and, of course, to obtain a connexion formula to 
link the two forms of solution. After several vain attempts, 
he returned to the problem in 1857, and after three days of 
concentrated effort the light (metaphorically) dawned at 
3 o'clock in the morning. Excitedly he wrote to his fiancee 
on 19 March, telling her of his new realization, but pref- 
acing it with a poignant acknowledgement of his under- 
standing that once they were married he would no longer 
be permitted to work to such hours 1 (Larmor, 1907). 

The matter is well illustrated by studying Airy's equa- 
tion: 



dx 2 



— xy = 



(49) 



which has a turning point at x = 0, and whose solutions 
can be expressed in terms of Bessel functions of order |. 
The advantage of using such a simple equation, which, in- 
cidentally, is central to the discussion of the more general 
equation (1) which follows, is that there are exact integral 
representations of its solutions which are analytic across 
the turning point. These are 



l 

2ni 



dt 



(50) 



the integral being along an appropriate contour C in the 
complex-^ plane, which canonically is one of the curves 
C_, Co or C+ depicted in Fig. 1. There are various ways 
of deriving this result, but they are not material to my 



presentation here. All that one need do is to substitute 
the representation (50) into equation (49) and show that 
it fits. There arc three natural solutions, y_, j/o an d y+, 
represented by equation (50) on the contours C_, Co and 
C + , respectively. Of course they cannot be independent, 
because equation (49) is of only second order. Since the 
integrand is entire, it follows that the integral along C_ U 
Co UC+ vanishes, and therefore y-(x) + yo(x) + y+(x) = 
for all (finite) x. 

To span the solution space it is customary to adopt the 
two independent functions: 

Ai(x) = y (x) , (51) 

evidently named after Airy (by Harold Jeffreys), and 



Bi(x) = iy-(x) - iy+{x), 



(52) 

a not unnatural consequent appellation. Both functions 
are real for real x. It is Ai(x) which is of principal interest 
here, for that is the solution that is relevant for our main 
purpose, for it decays, exponentially, for large positive 
(real) x. It can be obtained in terms of a single integral 
J of a complex variable by deforming the contour Co to 
(-S+) U 5„, on which t = e'^^s or t = e+"/ 3 s with s 
real, whence 



yoix) = h 



e iir/3j( e iT/3 x j _ e -iir/3j( e -iT/3 x j 



which is evidently real when x is real; 
I{z): = 



(53) 



(54) 



The function Bi may be expressed similarly. 

Series expansions of the Airy integral (50) were devel- 
oped by Stokes (1864, 1871); Jeffreys and Jeffreys (1956) 
describe how to obtain asymptotic approximations for 
large \x\ by Debye's method of steepest descents. When 
x > 0, the saddle points in the complex- i plane are at 
t = ±\fx, only the one on the positive real axis being ac- 
cessible to Co. In its vicinity the line of steepest descents 
is parallel to the imaginary axis, and by deforming Co to 
pass through t = ^fx along that line one obtains immedi- 
ately for Ai(:r): 



ki{x) 



exp 



(-fxi)[l + 0(^)] (55) 



as x — * +oo. If x < 0, the saddles are at t = ±i\/—x, the 
lines of steepest descents being inclined at ±7r/4 from the 
real axis, respectively. Neither is accessible to Co, but we 
may evaluate y± instead and express yo as minus their 
: (^V 1 (-x)-V*exp [±i (|(-x) 3 / 2 + I)] , 

)^^-5(_ x )-i siri (|(_ a; )l + i 7r )[i + 0(x- 1 )](56) 



sum: y± 
whence 



Ai(x 



as x — > —oo. 

The solution Bi(x) grows exponentially as x — » oo, as of 
course does any other combination of j/_ , yo and y+ that 
contains a nontrivial component of either j/_ or y+. The 
particular combination (52) is chosen for defining Bi(z) 



except occasionally 
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because it contains no exponentially small component in 
its asymptotic expansion as x — ► +00. 

The asymptotic analysis of Bi(x) for x > requires anal- 
ysis in the vicinity of the saddle at t = —\/x, through which 
the line of steepest descents is along the (negative) real 
axis. Both C + and C_ can be distorted to pass through it, 
which leads to a doubling of the amplitude factor: 



Bi(z) 



' 7r 2 x 4 exp 



3 

X2 



[l + OOr- 1 )] 



(57) 



as x — > +00. For x < one appropriately combines the 
expressions for y_ and y+ obtained previously to yield 

Bi(z) - n~ 3 (-x) - 3 cos (| (-a;) * + ±tt) [1 + O (aT 1 )] (58) 

as a; — > —00. Notice that the Wronskian of the asymptotic 
representations (55) and (57), and (56) and (58) takes the 
same constant value, 7r _1 , either side of the turning point, 
as it should. 

One can develop expansions to higher order, but I do 
not do so here. I simply point out that for sufficiently 
large x the relatively small correction to expression (57) 
at any order exceeds even the leading-order expression 
(55). Stokes' realization was that it is because of that 
that one cannot analytically continue the leading-order 
expressions around the turning point and expect them still 
to represent the same solution of the differential equation 
(49). 

Note that expression (56) provides the appropriate 
phase of the oscillatory branch of the solution in this rel- 
atively simple case, which was Stokes's goal. But it does 
not yet answer the question for the more general equa- 
tion (1). However, one should perhaps pause for a mo- 
ment to appreciate the power of the argument in the Airy 
case. Although the asymptotic expressions (55)-(58) are 
necessarily only approximations, they are approximations 
to the exact representation (50), and the connexion be- 
tween them is therefore robust. It isn't even necessary to 
know what the functions look like for moderate or small 
values of x, where the conditions motivating the asymp- 
totic expansions at large \x\ are not satisfied, although 
as a matter of natural curiosity one might wish to know. 
Accurate numerical solutions satisfy that desire. With 
such a connexion securely understood, Rayleigh (1912) 
applied the analysis of the Airy integral to a study of the 
reflection of waves propagating through a medium (actu- 
ally a stretched membrane, for which equation (7) holds 
with u> c = 0) in which c -2 varies linearly with z. 

5.2. Jeffreys' connexion 

The connexion between the asymptotic representations 
of the solutions to the more general equation (1) either 
side of a turning point was first established by Jeffreys 
in his Cambridge Adams prize essay in 1923 (see also 
Jeffreys, 1925). Subsequently, interest in the issue arose 
in quantum mechanics with respect to solutions of the 
time-independent Schrodingcr equation with a Coulomb 
potential, in studying, in particular, the asymptotic en- 
ergy levels of the hydrogen atom. Brillouin (1926) re- 



lated Schrodingcr's equation to Hamilton- Jacobi theory 
of classical mechanics, and Wentzcl (1926) discussed the 
leading-order LG expansion of Schrodingcr's equation; in 
response, Kramers (1926), who as an impoverished young 
man lived as a guest of the Jeffreys in their house in 
order to make ends meet during a stay in Cambridge, 
added the connexion formula that Harold Jeffreys had es- 
tablished. Mathematically their analysis was not new, 
but, being applied to a branch of physics that was more 
fashionable than the classical wave propagation that had 
interested Rayleigh and Jeffreys, it attracted the atten- 
tion of more physicists (e.g. Young and Uhlcnbcck, 1930; 
Kemble, 1935), and also mathematicians such as Langer 
(1934, 1937, 1949), and appellations combining the ini- 
tials of the surnames of the three quantum physicists, in 
various orders, were given to the method, eventually con- 
verging on the order WKB. Yet later, when Jeffreys' pi- 
oneering work was more widely recognized, the initial J 
was added, either in front or, more commonly, behind. I 
adopt the former, partly because Jeffreys has precedence, 
and partly because his analysis was designed to solve a 
wider class of problems than merely those associated 
with the time-independent 

Schrodinger equation. It is interesting to contemplate, 
however, that had those who originally named this om- 
nipresent and forgiving approximation been more aware 
of its true history, they might have dropped the initials 
W, K and B, and instead called it the SJ approximation 2 . 

There are now several ways of justifying the approxi- 
mation (e.g. Heading, 1962). Jeffreys (1925) noted that 
if f(z) has a simple zero at z = Zq, namely /(zo) = and 
f (z ) ^ 0, then f(z) ~ f'(z )(z - z ) near z = z ; and 
if f'(zo) > 0, the substitution x = —k^[f'(z ))3(z — z ) 
transforms equation (1) approximately into Airy's equa- 
tion (49), which enables approximate solutions of the full 
equation to be likened to the exact solutions of the com- 
parable approximate equation, and hence to a connexion 
between representations of the solutions either side of the 
turning point (via the integral representation (50)). He 
pointed out that the outcome is a limiting form of the 
true solution if k is arbitrarily large. Thus the phase of 
that oscillatory solution in z < z that connects to the so- 
lution that decays as z — > +00 is determined from equation 
(56). 

Nowadays it is customary to use an approximation that 
is also valid far from the turning point. This is achieved by 
applying what has become known as the Liouville-Green 
transformation to equation (1), namely 



x = -fc 3 sgn(/) 



p-dz 



* = (7$ := {—xf 



(59) 
(60) 



which leaves the equation in normal form. It may be writ- 
ten as 



4-^- - x$ = k 2 /i(x)$ , 
ax z 

after Stokes and Jeffreys 



(61) 
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where 

h{x[a(z)}} = -a 3 ^. (62) 

First it should be noted again that near the turning point 
/ ~ /q(z — zo) , where /q = f'(zo), and that therefore 
x ~ — sgn(/Q)fc3 3 (z — zo), which is Jeffreys' transfor- 
mation: the term on the right-hand side of equation (61) 
arises from small higher-order terms in the Taylor expan- 
sion of / about zo, and can be neglected. Far from the 
turning point / is no longer small, but, by definition, is of 
order unity. One can therefore estimate the ratio of the 
right-hand side of equation (61) to the second term on the 
left-hand side to be 0(ArV 4 /a;#|) = 0(fc- 2 i?/ 2 / -1 ) = 
0{k~ 2 HJ 2 ) = o(l) as k — > oo. Therefore, provided that 
k is sufficiently large, the right-hand side can be ignored 
throughout, reducing equation (61) to Airy's equation. 
What I have described here is not a proof, but it suggests 
that the solution, $ = Ai(x), to what is now called the 
comparison (Airy) equation provides a valid approxima- 
tion to the required solution of equation (1). And indeed 
Olver (e.g. 1974, 1978) has shown that for sufficiently 
well behaved functions / the Airy-function approximation 
converges uniformly for all x to the solution of equation 
(61) as k — ► co. 

It should be recalled that the JWKB procedure de- 
scribed above applies to simple turning points, for which 
f'(zo) + 0. If f (z ) = but f"(z ) ^ 0, then one can carry 
out a similar analysis using Weber's equation in the form 
d 2 y/dx 2 

+x 2 y = as the comparison, and one may generalize fur- 
ther to yet higher-order turning points. Olver has proved 
uniform convergence in such cases too. 

The uniformly valid JWKB approximation to the solu- 
tion to equation (1) having a single simple turning point 
at z = zo such that / > for z > zo can therefore be writ- 
ten in terms of Ai(— x), where x is defined by equation 
(59). Far from the turning point the appropriate asymp- 
totic approximation (55) or (56) to the Airy function is 
applicable. It yields, after setting / = k 2 , as in §1, and 
relating $ to iff by equation (60), 

i f z ° 

iff ~ iA|/t|~2 exp(- / \n\dz) for z < z , (63) 

J z 

and 

* — An~i sin (^J ndz+^j for z>z , (64) 

where A is a constant amplitude. These are equivalent to 
the Liouville- Green expressions (19) and (16); but they 
are more precise, because they define which combination 
of the solutions (16) corresponds to the evanescent solu- 
tion in z < zq. It goes without saying that there is a similar 
pair of JWKB-approximate solutions that corresponds to 
exponential growth in the so-called evanescent zone. Such 
solutions would need to be combined with (63) and (64) if 
one were solving, for example, a barrier-penetration prob- 
lem, either classical or quantum. 



6. Equations with two well separated turning 
points: eigenvalue problems 

Few sharp boundaries are encountered in astrophysics. 
Normally, in wave problems, when waves are confined 
within some region V it is because they encounter other, 
bounding, regions £ into which they cannot propagate, 
but the 'interface' is of finite thickness and is smooth. (I 
have in mind a 3-dimensional space with propagating and 
evanescent regions V and £ separated by 2-dimensional 
surfaces.) That is not to say that no wave at all can prop- 
agate in the bounding regions £. Usually the bounding 
regions can support waves of the same type as those un- 
der consideration, but with rather different frequency or 
wave number. The 'boundary' B between those regions 
can usually be regarded as the location of a turning point 
(or, more generally, as the locus of a turning point), in 
the sense that I have used it in this article, and is as well 
defined as arc the critical cutoff frequencies. 

I confine my discussion to situations in which the back- 
ground state is independent of time, so that frequency 
uj is well defined, and is conserved, as is wave energy 
in a dissipationless system. (An LG expansion in time, 
and even a JWKB approximation, can be made to study 
waves in an appropriately slowly changing environment, 
but I do not address that here.) Consider there to be 
a set of locally Cartesian coordinates (£,r?,C) established 
in the vicinity of the boundary, with £ perpendicular to 
the boundary. (The co-ordinate £ here is not to be con- 
fused with the vertical component of the displacement of 
equation (33).) Generically, the background (equilibrium) 
state varies more rapidly with £ than it does with £ and 
r\ (there are exceptions). On a length scale smaller than 
the scale of variation on B of the background state, the 
differential equation describing the waves reduces approx- 
imately to an ordinary (linear) differential equation with 
respect to C (often of second order, analogous to equation 
(1), the only case I consider explicitly here), with a turn- 
ing point 'on' B. Why does that turning point arise? In 
other words, why is wave propagation possible in V but 
not in £? 

One reason might be that no wave of frequency uj can 
propagate in £, whatever its orientation. Then a wave 
incident on B has its direction reversed - equation (1) 
permits no other possibility - undergoing what is called 
(true) reflection. Alternatively, £ simply doesn't permit 
continued propagation of a wave with a particular angle 
of incidence: if the background state is independent of 
some particular coordinate £ (or 77, it cannot be 0, then 
the component k^ of the local wavenumber is conserved 
across the boundary, and if / in equation (1) is then ren- 
dered negative in £ for that value of k^ (at the frequency 
uj of the wave), propagation is impossible; indeed, because 
the variation of the background state is gentler in B than 
perpendicular to it, that describes approximately the sit- 
uation with respect to the entire component k± of the 
wavenumber in B (i.e. perpendicular to the normal). The 
process is often called total internal reflection, and, if / is 
continuous, the case first considered by Rayleigh (1912), 
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can be thought of as the result of continuous refraction 
in the vicinity of B. Note that because the exponents of 
the two LG solutions (16) have the same magnitude at 
any point, the angle of reflection is equal to the angle of 
incidence. Note also that, because / varies when either ui 
or k± vary (at given £), the location of B varies with not 
only frequency but also with the angle of incidence of the 
wave; and, for some values of lo and fcj_, B may not exist. 

A region V that is enclosed by 8 is called a cavity. Any 
wave confined within it, in the absence of dissipation, will, 
under most situations, eventually pass arbitrarily close to 
any point it had passed formerly, and interfere with itself. 
For certain frequencies, that interference is constructive 
everywhere in V, and resonance occurs. The resulting 
motion is called a mode of oscillation: the contents of the 
entire region V oscillate in unison with frequency lo (it 
being quiescent elsewhere). For sufficiently high k the LG 
expansion can be used along the phase trajectories of the 
waves, using JWKB theory in the vicinity of the boundary 
of V, to determine the values of lo that permit resonance. 
These values are called the eigenfrcqucncies of the cavity. 
Such an analysis can be geometrically complicated enough 
to detract from the main point of this discussion, so I post- 
pone discussion of the general problem to another article. 
Here I simply consider a problem in just one dimension, 
in which a wave after reflection is bound to pass through 
all the points it had passed through previously. Then the 
governing differential equation is of the type (1), with k 
and / depending in some way on to. 

To be more specific, let us consider a situation in 
which f(z) = k 2 > for z\ < z < 22, and / < elsewhere. 
Asymptotic solutions ^2(2) in the vicinity of 22, and far 
from z±, may be written as a{x)<S>2{x), where x is defined 
by equation (59) with zq replaced by 22, <J is defined by 
equation (60), and 3>(x) = Ai(x). In the vicinity of z\ one 
may do likewise to define ^i(^) = cr(x)$>i(x), but with x 
defined with the opposite sign, and now with zq replaced 
by z\ , to ensure that the interval in which Ai is oscillatory 
satisfies z> Z\. For simplicity I shall assume that z\ and 
Z2 are sufficiently far apart that there is a common region 
of validity in (21,22) of the high-(— x) expansion (55) of 
<£>i and <f>2- This must be possible for sufficiently large k. 
Then "fi and ^2 can be matched. According to equation 
(64), 



which is possible only if 



2 s i n J k J ndz + jj for z » 21 

(and 2 < 2 2 ) (65) 



and 



^ 2 ^A 2 n 3 sin [k / red2 + f for z < z 2 



(and x>zi) . (66) 

The representations are identical throughout the common 
interval of validity if A2 = ±Ai and 



cos 



\k / Kd2+f I =0 



(67) 



(68) 



where n is an integer, and, in order to achieve a possible 
matching with the two evanescent solutions either side 
of the oscillatory region, at least one half-wavelength is 
required: therefore n = 1,2,3, ... . This is the desired 
eigenvalue equation. 

I promised to illustrate the approximation with 
acoustic-gravity modes of a star, for which, in the case 
of a spherical star, kn is K given by equation (37). I do 
so with some hesitation, however, because stars are not 
quite as simple as the situation I have just discussed. The 
reason is that if lo is high, which is one way of making 
K large, then the upper turning point r2, which occurs 



roughly where lo 2 (r) 



is close to the surface, towards 



which lo 2 , and perhaps also N (depending on whether or 
not the outer layers are convective, although the lo 2 di- 
viding iV 2 makes its effect relatively small), rise rapidly, 
as though they are approaching a singularity located at 
what I call the seismic surface of the star. This is a struc- 
tural property of all stars. In practice, where the outer 
layers become optically thin, the stratification becomes 
approximately isothermal, and therefore essentially expo- 
nential, and the singularity that would be encountered 
in a polytrope, for example, is avoided. Nevertheless, for 
many modes, with lo rather less than the value of tu c in the 
photosphere, the transition to exponential behaviour oc- 
curs sufficiently well beyond the upper turning point for 
its presence to be hardly discernible in the dynamically 
active propagating regions below. Therefore the solution 
to equation (1) actually feels the influence of the phan- 
tom singularity, and the eigenvalue equation (68) needs 
adjustment. How that is achieved is beyond the scope of 
this elementary introduction, and is postponed to a subse- 
quent article. Here I simply assume that lo is not so high 
as to make that adjustment necessary. There is also a sim- 
ilar problem near the centre, where L 2 /r 2 suffers a true 
(co-ordinate) singularity. This is a geometrical property of 
any wave problem in a sphere, and was first encountered in 
the context of JWKB theory in quantum mechanics (with 
a Coulomb potential). That singularity is well avoided 
if the degree I of the spherical harmonic (the angular- 
momentum quantum number in quantum mechanics) is 
sufficiently large, so that the lower turning point ri, which 
is determined approximately for high-frequency (acoustic) 
modes by ri/c(ri) = L/lo, is not very close to the centre. 
Similar problems face gravity modes, which have very low 
frequency, although the problem near the surface of the 
star is present only in stars with radiative envelopes, for 
otherwise the modes are shielded from the surface by the 
evanescent convection zone in which N 2 < 0. 

With these caveats in mind, I now intrepidly make the 
substitution kn = K in equation (68), using equation (37) 
to obtain 



U 2 



dr ' 



(69) 
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hoping that a range of uj exists which is high enough 
for the LG expansion to be reasonably accurate and low 
enough for the JWKB approximation near the turning 
points not to be unduly disturbed by singularities, phan- 
tom or real. 

For acoustic modes uj 2 is large, and for many purposes 
N 2 /uj 2 may be neglected compared with unity, yielding 



, co 2 



dr 



(70) 



where n = t{t{), = r(ra) and w = uj/L; for gravity 
modes uj 2 /c 2 may often be neglected, yielding 



(n-i)Tr 



N 



- 1 - 



L 2 c 



dr 



(71) 



in both cases the turning points r\ and r^ being inter- 
preted as the points at which the integrands vanish. 

It is interesting to develop equation (70) further in the 
case when L is large, for then the lower turning point is 
close enough to the surface for the effect of the spheri- 
cal geometry to be small. In that case one may regard 
the horizontal phase velocity v<f,h = rw to be essentially 
independent of r. If, in addition, one approximates the 
surface layers of a star of seismic radius R by a plane- 
parallel polytrope of index so that c 2 = jigz/ (p + 1) 
and uj 2 = jti(/i + 2)7ig , /4(^+ l)z, where z — R — r is depth 
beneath the seismic surface of the star (e.g. Gough, 1993), 
then the acoustic depth f beneath the seismic surface is 
given by f = 2z/c, and uj 2 = fi(fi + 2)/f 2 . Equation (70) 
reduces to 



fi(fi + 2) 



iWr 2 



4(/i- 



l) 2 v 2 



Mi 



df 



(M+%h y^7+2)-2 



2aJ 



(72) 



in which f\ and are respectively the acoustic depths of 
the upper and lower turning points. This equation may 
be rewritten 



(n + a)n 



= F(v 4>h ) , 



(73) 



where F(v) := (p + l)iro/2jig and a = \ [a/ fi(fi + 2) - 1] = 
constant, which is a special case of what is now known 
as Duvall's law (1982), namely equation (73) with the 
functional form of F unspecified and with the constant a 
not necessarily related directly to a polytropic index, and 
which formally describes high-frequency acoustic modes 
in a stellar envelope with a (well behaved) reflecting sur- 
face, whatever its stratification. The polytropic form de- 
rived here can be rewritten ui 2 ~ (n + a)uj^Rk as k — > oo, 
where ujq = 2jig/(n+ 1)R is the square of a characteris- 
tic acoustic frequency of the outer layers of the star, and 
here k = L/R is the horizontal wavenumber at the sur- 
face. This is the parabolic form which approximates the 
solar k — oj relation familiar to hclioseismologists, and was 
used to provide the first seismic calibration of the adi- 
abatic constant deep in the solar convection zone. The 



more general form of equation (73) when the polytropic 
sound-speed variation is not assumed can be inverted to 
give c/r as a function of the observables [n + a)ir/uj and 
w (e.g. Gough, 1993), and was used to provide the first 
inferences of the sound speed through the Sun. 

7. Closing words 

The JWKB approximation provides a wonderfully ro- 
bust representation of the waves encountered throughout 
physics. Turning points are very common, and, provided 
one is well clear of them, the approximation reduces to 
simple formulae in terms of an exponential function on 
one side, and a trigonometrical function on the other, to- 
gether with the connexion formula providing the value 
of the phase of the trigonometrical function according to 
whether the exponential function grows or decays away 
from the turning point. Even near the turning point, it is 
straightforward to evaluate the Airy-function representa- 
tion from which these simpler formulae are derived. 

In solar physics, the leading term in the LG expan- 
sion, now commonly called the WKB approximation, is 
almost ubiquitous in studies of waves in the atmosphere, 
and when turning points are encountered the more pow- 
erful JWKB approximation is available. For internal 
acoustic-gravity waves at least one turning point, r\, is 
always present, and the connexion formula provided by 
the JWKB approximation is essential. I have discussed 
the latter explicitly, showing how, if uj is low enough for 
the waves to be trapped inside the Sun by reflection near 
the surface, the approximate equation (69) determines the 
eigenfrcquencies of the resonant modes. In the case of 
acoustic modes of the Sun and Sun-like stars, that formula 
is only approximate. There are two reasons: (i) if uj is low 
enough for the upper turning point not to be too close to 
the phantom singularity at the seismic surface of the Sun, 
then one runs the risk of not having k large enough for 
the Airy equation (61) with the right-hand side ignored to 
provide a reliable comparison to the full equation. That 
risk is usually small, because in most (but not all) cases 
the Airy- function representation is amazingly robust; (ii) 
if the frequency is high, not only does the upper turn- 
ing point approach the acoustic surface of the Sun, but it 
also encounters the upper superadiabatic boundary layer 
of the convection zone where, over a small distance, the 
acoustic cutoff frequency becomes imaginary; then there 
appear to be three turning points which arc too close to- 
gether for the technique described in §6 for piecing to- 
gether different Airy functions to be applicable. In that 
case special attention beyond the scope of this introduc- 
tion is required, but I might at least point out that it is 
pertinent to the interpretation of observational investiga- 
tions of the transition to supercriticality, namely to the 
high-frequency domain exceeding the acoustic cutoff fre- 
quency in the atmosphere. In the case of gravity waves in 
a star with a convective envelope, the threat of a singular- 
ity occurs only near the lower turning point, although in 
the absence of a convective envelope, the surface threatens 
too. 
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It is worth mentioning that the leading-order LG expan- 
sion is formally applicable only in the limit of infinitesi- 
mally slow variation of the function / in equation (1). It 
therefore does not capture reflection. The JWKB approx- 
imation behaves likewise, except near the turning point. 
Reflection in a spatially varying propagating region is 
likely to occur everywhere to some degree, but it is signifi- 
cant only when the scale of variation Hf of the background 
state is comparable with or much smaller than the charac- 
teristic inverse wavenumber (kf? )~ 1 of the wave. When it 
is much smaller the variation may be regarded as a discon- 
tinuity, across which one can match two JWKB solutions, 
in the manner adopted by Poisson (1817) under condi- 
tions when the simpler approximation obtained by taking 
/ to be piecewise constant is applicable. If kfzHf « 1 
over an extended region, one might need to adopt a more 
complicated comparison equation for the JWKB approxi- 
mation: to deal with a simple smooth barrier between z\ 
and 22) for example, one can adopt Weber's equation in 
the form y" + k 2 (z~ z\){z — z 2 )y = (e.g. Langer, 1959); 
more complicated variations in / require correspondingly 
more complicated comparison equations. One must con- 
sider in such cases whether the result would justify the 
effort, because direct numerical solution of an ordinary 
differential equation is much more straightforward, even 
when there are singularities present. 

Whether the result justifies the effort depends on the 
use to which one wishes to put the result. If all one needs 
are eigenfunctions and their corresponding eigenvalues, 
the direct numerical computation is in many cases simpler 
and more reliable. But analytical results are, for many a 
person, easier to interpret. And then it is more likely that 
one could be able to put them to good use. For exam- 
ple, in the early days of hclioscismology, the significance 
of the so-called large and small frequency separations of 
low-degree acoustic modes was first recognized from an 
asymptotic analysis, as subsequently, in a sense, was the 
Duvall law (1982), use of which has been central to many 
helioseismological investigations. It has been argued by 
some helioseismologists that asymptotic analysis was ac- 
tually unnecessary because the discoveries could equally 
have been made by numerical investigation. Perhaps they 
could. But the fact of the matter is that they were not, 
and that progress in the subject would undoubtedly have 
been slower without asymptotics. That mode of discov- 
ery is likely to continue: the days of asymptotic analysis 
in stellar physics are certainly not gone. 

Acknowledgments 

I thank J. B. Keller for useful discussion, and 
J. Christenscn-Dalsgaard for the hospitality of the 
Department of Physics and Astronomy, University of 
Aarhus , where this article was written. Support from the 
European Helio- and Asteroseismology Network (HELAS) 
for participation at the HELAS workshop in Nice is grate- 
fully acknowledged. HELAS is funded by the European 
Union's Sixth Framework Programme. 



References 

Booker, J.R., Bretherton, F.P.: 1967, J. Fluid Mech. 27, 513- 
539 

Brillouin, L.: 1926, Comptes Rend. Acad. Sci. Paris 183, 24- 
26 

Carlini, F.: 1817, Giorn. Fis. Chim. Storia Nat. Med. Arti 10, 
458-460 

Christensen-Dalsgaard, J., Cooper, A.J., Gough, D.O.: 1983, 

Mon. Not. Roy. Astron. Soc. 203, 165-179 
Deubner, F-L., Gough, D.O.: 1984, Ann. Rev. Astron. 

Astrophys. 22, 593-619 
Duvall Jr, T.L.: 1982, Nature 300, 242-243 
Froman, N., Froman, P.O.: 1996, Phase-Integral Methods 

Allowing Near-lying Transition Points, Springer Tracts in 

Nat. Phil. 40, Springer, Heidelberg 
Gough, D.O.: 1993, in Astrophysical Fluid Dynamics (ed. J-P. 

Zahn & J. Zinn- Justin) , Elsevier, Amsterdam 399-560 
Green, G.: 1837, Trans. Camb. Phil. Soc. 6, 457-462 
Heading, J.: 1962, An Introduction to phase-integral methods, 

Methuen, London 
Jeffreys, R: 1925, Proc. Lond. Math. Soc. 2nd Ser. 23, 428- 

436 

Jeffreys, R, Swirles, B. (Lady Jeffreys): 1956, Methods 
of Mathematical Physics, Cambridge University Press, 
Cambridge 

Kemble, E.C.: 1935, Phys. Rev. 48, 549-561 

Kramers, H.A.: 1926, Zs Phys. 39, 828-840 

Lamb, R: 1909, Proc. Lond. Math. Soc. 2nd Ser. 7, 122-141 

Langer, R.E.: 1934, Bull. Am. Math. Soc. 40, 545-582 

Langer, R.E.: 1937, Phys. Rev. 51, 669-676 

Langer, R.E.: 1949, Trans. Amer. Math. Soc. 67, 461-490 

Langer, R.E.: 1959, Trans. Amer. Math. Soc. 90, 113-142 

Larmor, J. (arr.): 1907, Memoir and Correspondence by the 

late Sir George Gabriel Stokes, Bart, Cambridge University 

Press, Cambridge, p. 62 
Liouville, J.: 1837, J. Math, pures appl. 2, 16-35 
Olver, F.W.J.: 1974, Asymptotics and Special Functions, 

Academic Press, New York 
Olver, F.W.J. : 1978, Phil. Trans. Roy. Soc. London 289, 501- 

548 

Poisson, M.: 1817, Mem. Acad. Roy. Soc. Inst. France 2, 305- 
402 

Rayleigh: 1912, Proc. Roy. Soc. London A86, 207-226 

Reye, T.: 1872, Die Wirbelsturme, Tornados und 

Wettersaulen, Genesius, Halle 
Schmitz, F., Fleck, B.: 1998, Astron. Astrophys. 337, 487-494 
Stokes, G.G.: 1864, Trans. Camb. Phil. Soc. 10, 105-128 
Stokes, G.G.: 1871, Trans. Camb. Phil. Soc. 11, 412-425 
Wentzel, G.: 1926, Zs Phys. 38, 518-29 

Young, LA., Uhlenbeck, G.E.: 1930, Phys. Rev. 36, 1154- 
1167 



